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Abstract 

In this work, we propose a mixed finite element method for solving elliptic mul¬ 
tiscale problems based on a localized orthogonal decomposition (LOD) of Raviart- 
Thornas finite element spaces. It requires to solve local problems in small patches 
around the elements of a coarse grid. These computations can be perfectly par¬ 
allelized and are cheap to perform. Using the results of these patch problems, we 
construct a low dimensional multiscale mixed finite element space with very high ap¬ 
proximation properties. This space can be used for solving the original saddle point 
problem in an efficient way. We prove convergence of our approach, independent 
of structural assumptions or scale separation. Finally, we demonstrate the applica¬ 
bility of our method by presenting a variety of numerical experiments, including a 
comparison with an MsFEM approach. 

Keywords mixed finite elements, multiscale, numerical homogenization, Raviart-Tho- 
mas spaces, upscaling 
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1 Introduction 

In this work we study the mixed formulation of Poisson’s equation with a multiscale 
diffusion coefficient, i.e. where the diffusion coefficient is highly varying on a continuum 
of different scales. For such coefficients, the solution is typically also highly varying and 
standard Galerkin methods fail to converge to the correct solution, unless the features on 
the finest scale are resolved by the underlying computational mesh. A classical application 
is the flow in a porous medium, modeled by Darcy’s law. In this case, the multiscale 
coefficient describes a permeability field, which is heterogeneous, rapidly varying and 
has high contrast. Classical discretizations that involve the full fine scale often lead 
to a vast number of degrees of freedom, which limits the performance and feasibility 
of corresponding computations. In this paper, we address this kind of problems in the 
context of mixed finite elements. 

We will interpret the mixed formulation of Poisson’s equation in a Darcy flow set¬ 
ting, referring to the vector component as flux, and the scalar component as pressure. 
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In Darcy flow applications the flux solution is of particular interest since it tells us how 
a fluid is transported through the medium. It is desirable and common to use flux con¬ 
servative discretization schemes. The proposed method is based on the Raviart-Thomas 
finite element |3Tj which is locally flux conservative. Concerning the mixed formulation 
of Poisson’s equation, corresponding multiscale methods were for instance proposed in 
[hme]. These methods are based on the Raviart-Thomas finite element and fit into 
the framework of the Multiscale Finite Element Method (MsFEM, cf. [19]). Another 
family of multiscale methods is derived from the framework of the Variational Multiscale 
Method (VMS) [2U[ [2T, 22HMH25]- Multiscale methods for mixed hnite elements based on 
VMS are proposed and studied in [U ES, 2B] • Inspired by the results presented in [28], a 
new multiscale framework arose [25] • We refer to this framework as Localized Orthogonal 
Decomposition (LOD). It is based on the idea that a hnite element space is decomposed 
into a low dimensional space that incorporates multiscale features and a high dimensional 
remainder space which is given as the kernel of an interpolation or quasi-interpolation 
operator. The multiscale space can be used for Galerkin-approximations and allows for 
cheap computations. Various realizations have been proposed so far. For corresponding 
formulations and rigorous convergence results for elliptic multiscale problems, we refer 
to la m da na m\ for Galerkin hnite element methods, to [[131 ITT] for discontinuous 
Galerkin methods and to in for Galerkin Partition of Unity methods. Among the var¬ 
ious applications we refer to the realizations for eigenvalue problems [27], for semilinear 
equations [16], for the wave equation [3] and for the Helmholtz equation [30] . 

In this paper we introduce a two level discretization of the mixed problem, that is we 
work with two meshes: A fine mesh (mesh size h) which resolves all the fine scale features 
in the solution and a coarse mesh (mesh size H) which is of computationally feasible size. 
This gives us a fine and a coarse Raviart-Thomas function space for the flux. We denote 
them respectively by V), (high dimensional) and Vh (low dimensional). The kernel of the 
(standard) nodal Raviart-Thomas interpolation operator n# onto Vh is the detail space 
V[. This space can be interpreted as all fine scale features that can not be captured 
in the coarse space Vh- A low dimensional ideal multiscale space is constructed as the 
orthogonal complement to the divergence free fluxes in V We prove that this space has 
good approximation properties in the sense that the energy norm of the error converges 
with H without pre-asymptotic effects due to the multiscale features. However, the basis 
functions of the ideal multiscale space have global support and are expensive to compute. 
We show exponential decay of these basis functions allowing them to be truncated to 
localized patches with a preserved order of accuracy for the convergence. The resulting 
space is called the localized multiscale space. The problems that are associated with the 
localized basis functions have a small number of degrees of freedom and can be solved in 
parallel with reduced computational cost and memory requirement. Once computed, the 
low dimensional localized multiscale space can be reused in a nonlinear or time iterative 
scheme. 

We prove inf-sup stability and a priori error estimates (of linear order in H) for both 
the ideal and the localized method. The local L 2 -instability of the nodal Raviart-Thomas 
interpolation operator leads to instabilities as h decreases for the localized method. We 
show that these instabilities can be compensated by increasing the patch size or using 
Clement-type interpolators instead. In the numerical examples we verify that the local¬ 
ized method has the theoretically derived order of accuracy. We confirm our theoretical 
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findings by performing experiments on the unit square and an L-shaped domain, as well 
as using a diffusion coefficient with high contrast noise and channel structures. The pro¬ 
posed method is also compared numerically with results from an MsFEM-based approach 
using a permeability field from the SPE10 benchmark problem. 


2 Preliminaries 

We consider a bounded Lipschitz domain C (dimension d — 2 or 3) with a piecewise 
polygonal boundary <912 and let n denote the outgoing normal vector of <9f 2. For any 
subdomain w C we shall use standard notation for Lebesgue and Sobolev spaces, 
i.e. for r G [1, oo], L r (oo ) consists of measurable functions with bounded //-norm and 
the space H l (oo ) consists of L 2 -bounded weakly differentiable functions with L 2 -bounded 
partial derivatives. The full norm on H 1 (uj ) shall be denoted by || • whereas the 

semi-norm is denoted by | • | h 1 ^) '■= || V • 

For scalar functions p and q we denote by (p, q) u := j'^pq the L 2 -scalar product on oo. 
When u = 12, we omit the subscript, i.e. (p, q) := (p, q )q. For d -dimensional vector valued 
functions u and v, we define (u, v) w := f u- v with (u, v) = (u, v)q. Observe that we use 
the same notation for norms and scalar products in L 2 without distinguishing between 
scalar and vector valued functions. This is purely for simplicity, since the appropriate 
definition is always clear from the context. We use, however, bold face letters for vector 
valued quantities. 

In the following, we define the Sobolev space of functions with L 2 -bounded weak 
divergence by Ff(div, oo) := {v G [L 2 (oo)] d : V • v G L 2 (oo)}. We equip this space with 
the usual norm || • ||tf(div,«), where ||v||^ (diVjU)) := ||V • v|| 2 2(w) + ||v|| 2 2(w) . Additionally, 
for oo = 12, we introduce the subspace fh 0 (div, hi) := {v G H( div, hi) : v • n|gn = 0} of 
functions with zero flux on the boundary, where v-n|gQ should be interpreted in the sense 
of traces. We denote by L 2 (f2)/M := {q G L 2 (Q) : f n q = 0} the quotient space of L 2 (Q) 
by M. The continuous dual space of a Banach space X is denoted by X'. 


2.1 Continuous problem 


With these definitions we are ready to state the continuous problem, which is Poisson’s 
equation in mixed form with Neumann boundary conditions on the full boundary. 

Definition 1 (Continuous problem). Find u G V := H 0 (dw,Q), p G Q := L 2 (Q)/R such 
that 


(A 1 u,v) + (V • v, p) = 0, 

(V • u ,q) = —(/, q), 


( 1 ) 


for all v EV, q G Q. 

We pose the following assumptions on the coefficient and data. 

Assumption A (Assumptions on coefficients, data and domain). 

(Al) A G [L 00 (f2)] d>< ' i is a diffusion coefficient, possibly with rapid fine scale variations. 
Its value is an almost everywhere symmetric matrix and bounded in the sense that 
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there exist real numbers a and /3 such that for almost every x and any v G M d /{0} 


0 < a < 


(A(x) x v) ■ v 

V ■ V 


< f3 < oo. 


(A2) / G L 2 (Q) is a source function that fulfills the compatibility condition J Q f — 0. 

(A3) The domain Q is a bounded Lipschitz domain with polygonal (or polyhedral) bound¬ 
ary. 

We introduce the following bilinear forms and norms. Let 


a(u, v) := (A x u, v) and 6(v, q) (V • v, q) 


and, further, 


IM|y : = IM|tf ( div,n) and ||<?IIq := IMUw 

The energy norm is defined as the following weighted flux L 2 -norm, 

ll|v||| J := ||A“ 1/2 v||| 2(n) = a(v,v) 


• • • 111 1112 
The energy norm can be subscripted with a subdomain u> C hi, for example |||-||| , to 

indicate that the integral is taken only over that subdomain. 

The following lemma gives the conditions for existence and uniqueness of a solution 

to the mixed formulation in (|TJ) for subspaces V C V and Q C Q. This lemma is useful 

for establishing existence and uniqueness for all discretizations presented in this paper, 

since all presented discretizations are conforming. 


Lemma 2 (Existence and uniqueness of solution to mixed formulation). Let V C V and 
Q C Q. Denote by /C = {v e V : b(v,q ) = 0 Vg G <2}. If a(-, •) is coercive on K. with 
constant a > 0, i.e. a(v, v) > a||v||y for v G /C, and bounded with constant (3 > 0, i.e. 
|a(v,w)| < /3||v|| vll w llv f or v > w ^ V, and additionally &(•,•) is inf-sup stable with 
constant 7 > 0, i.e. 


inf sup 

q&Q vgV 


v llv||?|lo 


> 7 , 


then the problem a(u, v) + b(v,p ) — b(u,q) = ( f,q ) for all (v, q) G V x Q has a unique 
solution (u,p) G V x Q bounded by 


||u|| y < %j^\\f\\ L *m and ||p||q < 4ll/IU 2 (n)- 

Proof. See e.g. 0, Theorem 4.2.3]. □ 

Under Assumptions (Al)-(A3), the conditions for Lemma|2]are satisfied for V = V and 
Q = Q with a = a, (3 = (3 and 7 being a constant that depends only on the computational 
domain. The lemma then yields a unique solution to the continuous problem ([Tj). 
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2.2 Discretization with the Raviart—Thomas element 

Regarding the discretization, we introduce two conforming families of simplicial (i.e. tri¬ 
angular or tetrahedral) meshes (7 h} and {Th} of where h and H are the maximum 
element diameters. Throughout the paper we refer to Th as the fine mesh and to Th as 
the coarse mesh. Hence, we indirectly assume h <C H. We pose the following assumptions 
on the meshes. 

Assumption B (Assumptions on meshes). 

(Bl) The fine mesh Th is the result of one or more conforming (but possibly non-uniform) 
refinements of the coarse mesh Th such that Th D Th = 0. 

(B 2 ) Both meshes Th and Th are shape regular. In particular the positive shape regularity 
constant p for the coarse mesh Th will be referred to below and is defined as p = 
miiiTeTH d dtam^ where Bt is the largest ball contained in the element T G Th- 

(B3) The coarse family of meshes {Th} is quasi-uniform, whereas {Th} could be obtained 
from an arbitrary adaptive refinement. 

Remark 3 (Quadrilateral or hexahedral elements). Affine quadrilateral (or hexahedral) 
elements can also be used. However, the definition of the Raviart-Thomas element pre¬ 
sented below in this paper is based on triangular (or tetrahedral) meshes. 

We denote by t and T an element of T or Th, respectively. Similarly e and E denote 
an edge (for d = 2 ) or a face (for d — 3) of the elements of Th and Th- Further, n e 
(respectively n E ) is the outward normal vector of an edge (or face) e (respectively E ). 
We continue this section by discussing finite element discretizations using the two meshes. 

We denote all polynomials of degree < k on a subdomain u by P k (ui) and a d- 
dimensional vector of such polynomials by [P A ’(o;)] d . We introduce the // 0 (div, Q)-conform- 
ing lowest (zeroth) order Raviart-Thomas finite element. For each fine element t G Th 
and coarse element T G Th, the spaces of Raviart-Thomas shape functions are given by 

ITT h (t) = {v| t = [P °(t)} d + xF °(t)} and 
IZT h (T) = {v| T = [P°(T)] d + xP°(T)}, 

respectively, where x — (x\,... ,xfi) is the space coordinate vector. The Raviart-Thomas 
finite element spaces on Th and Th are then defined as 

v h = {v G i7 0 (div, Q) : v| t G IZT h {t) Vi G %} and 
V H = {v G R 0 (div, n) : v| T G 7 ZT H (T) VT G T H }- 

The degrees of freedom (in the coarse and fine Raviart-Thomas spaces) are given by the 
averages of the normal fluxes over the edges (respectively faces for d — 3). We denote the 
degrees of freedom by 

A e (v) := 7*7 / v • n e and N E (v ) := / v • n E 

\ e \ Je \ E \ JE 

for the fine and coarse discretization, respectively. The direction of the normal n e (respec¬ 
tively n^) can be fixed arbitrarily for each edge (respectively face). Here, N e and N E are 
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bounded linear functionals on the space W := i7 0 (div, Q) D L S (Q), for some s > 2 . Note, 
that the additional regularity (i.e. L S (Q) for s > 2) is necessary for the edge integrals to be 
well-defined (cf. [ 8 j). We introduce the (standard) nodal Raviart-Thomas interpolation 
operators IR : W —* 14 and U H : W —> Vh by hxing the degrees of freedom in the natural 
way, i.e. II h and 11 ^ are defined such that 

ARIRv) = N e (v) and N E (U H v) = N E (y). 

Additionally, we let Qh C Qh C Q be the space of all piecewise constant functions on 74 
and Th with zero mean. We denote by Ph and Ph the L 2 -projections onto Qh and Qh, 
respectively. Using the fine spaces, we define the fine scale discretization of ([Tj) , which 
will be referred to as the reference problem. 

Definition 4 (Reference problem). Find U/, G 14 and ph G Qh, such that 

a{u h ,v h ) + b(v h ,ph) = 0 , 

b(u h ,q h ) =-(f,q h ), 


for all v h G 14 and q h G Q h . 

A similar problem can be stated with the coarse spaces VA and Qh with flux solution 
u h- The remainder of this section treats only the fine discretization. However, all results 
hold also for the coarse discretization. 

We denote the space of divergence free functions on the fine grid by 


K h := {v G 14 : V ■ v = 0}. 


(3) 


Remark 5 (Kernel of divergence operator). A natural definition of Kh for our purposes 
is Kh = {v G 14 : (V • v, qQ = 0 VR G Qh}- However, since we have V • v G Qh for all 
v G 14 (due to the definition of the Raviart-Thomas element), we can characterize Kh 
equivalently as done in (|3]). 

To establish existence and uniqueness of a solution to the reference problem, we use 
that R is divergence compatible, i.e. we have the commuting property V • n^v = Ph V ■ v 
for v G W, and that is bounded on W (but not on 14), i.e. there exists a generic 
/i-independent constant CV such that ||n^v||y < CRIMR f° r v £ W. Using this, the 
inf-sup stability of &(•, •) with respect to 14 and Qh follows: For q G Qh, 


sup 

veu h 


b(v,q) 


\v 


sup 

vew 


(V • n fe v, q) 
||n,v|R 


. (V- 

> sup 771 
vew L'W 


v,g) 


I w 


(V ■ w ,q) > (q,q) 

C'w||w||w — CwCn\\q\\ L 2(ty 


Cw'Cn'hWvm, 


(4) 


where w G W is chosen such that V ■ w — q and ||w|R < CAIAHziq). This is possible by 
letting w = V0 for a solution 0 to A0 = q in O with homogeneous Neumann boundary 
conditions. Now, applying Lemma [2] with V = 14, Q = Qh, /C = Kh, we can derive the 
constants a = a, fi = (3 and 7 = 7 := and establish existence and uniqueness 

of a solution to the reference problem ([2]). Note that the inf-sup stability constant 7 is 
independent of h and hence also holds for the pair of spaces Vh and Qh- 
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In the following, we are mainly interested in approximating the flux component U/,, 
of the solution. We treat as a reliable reference to the exact solution. Note that the 
L 2 -norm of the divergence error is controlled by the data 


V • u - V • u h || L2(Q) < ||/ - P h f || L 2 ( 0 )- 


For the energy norm of the flux error, we have the following error estimate in the energy 
norm for the lowest order Raviart-Thomas element: 

|||u - Uft.||| < C7i|u|tfi ( fi), 

where C is independent of h. For a problem with a coefficient A that has fast variations 
at a scale of size e, we have in general that |u|#i(Q) ~ e -1 . Hence, we require h • C e before 
we can observe the linear convergence in h numerically. We call the regime with h > e 
a pre-asymptotic regime. The goal of this work is the construction of a discrete space 
which does not suffer from such pre-asymptotic effects triggered by A. In the following, 
we assume that the foie mesh is foie enough (in the sense that h <C e) so that |||u — u^||| is 
sufficiently small and hence U/, a sufficiently accurate reference solution. With the same 
argument, the accuracy of the coarse solution u h will not be satisfying as long as H > e. 
Note that reference problem (|2j) never needs to be solved. It just serves as a reference. 

In the next section, we will construct the ideal multiscale space of the same (low) 
dimension as Vh, but which yields approximations that are of similar accuracy as the 
reference solution u/j (in particular in the regime H e). Throughout the paper, we do 
not consider errors that arise from numerical quadrature. For simplicity, we assume that 
all integrals can be computed exactly. 

3 Ideal multiscale problem 

In this section, we construct a low dimensional space that can capture the fine scale 
features of the true multiscale solution. We focus on constructing a good multiscale 
representation of the flux solution u only. We call it ideal since the reference flux solution 
is in this space for all / G Qh • This should be contrasted to a localized multiscale space 
to be introduced in Section [4j In addition to the spaces 14 and Vh defined above we 
introduce the following detail space as the intersection of the fine space and the kernel of 
the coarse Raviart-Thomas interpolation operator, 

bf = {v G 14 : n^v = 0}. 

Since Vjf is the kernel of a projection, it induces the splitting 14 = Vjy © V£, where Vh is 
low dimensional and V/ is high dimensional. We refer to V/ as the detail space. In this 
section we aim at constructing a modified splitting, where Vh is replaced by a multiscale 
space which incorporates fine scale features. 

3.1 Ideal multiscale space 

We will construct the ideal multiscale space by applying fine scale correctors to all coarse 
functions in Vh, he. so that (Id — G/)(I4d is the desired multiscale space for a linear 
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corrector operator Gh■ The corrector operator is constructed using information from the 
coefficient A, and has divergence free range in order to keep the flux conservation property 
of the coarse space. 

The definition of the corrector requires us to construct the splitting K h = K H © K f h 
with 

K{ := {v G K h : fl^v = 0}, and K H := Range((Il//)|tfJ. 

Next, we introduce an ideal corrector operator. We distinguish between local (element¬ 
wise) correctors and a global corrector. 

Definition 6 (Ideal corrector operators). Let a T (u, v) := (A _1 u, v)t for T G Th■ For 
each such T G 74, we define an ideal element corrector operator Gf : V —* K[ by the 
equation 

a(G^v,v f ) = a T (v,v f ) (5) 

for all v f G K { h . Furthermore, we define the ideal global corrector operator by summing 
the local contributions, i.e. Gh J2 tgt h ^h- 

The ideal corrector operators are well-defined since equation ([5]) is guaranteed a unique 
solution by the Lax-Milgram theorem due to the coercivity and boundedness of a(-, •) 
on Kfi Using the ideal global corrector operator, we can define the discrete multiscale 
function space by 

V™:= (Id-G h )(V H ), 

where Id is the identity operator. This space has the same dimension as Vh- Furthermore, 
it allows for the splitting 14 = Vjfh © F f. Note that the ideal multiscale space is the 
orthogonal complement of K f h with respect to a(-, •), i.e. 

«(v^W f ) = 0 (6) 

for all G Vfif and v f G K f h . 


3.2 Ideal multiscale problem formulation 

In this section, we use the previously defined ideal multiscale space to define a (prelimi¬ 
nary) multiscale approximation. The ideal multiscale problem reads as follows. 

Definition 7 (Ideal multiscale problem). Find G Vffi s h and pn € Qh, such that 


a(u-H% v h) + b(y h ,p H ) = 0, 

H u h%(1h) = — (/> Qh), 

for all v/j G Vffi s h and qn G Qh- 

Lemma 8 (Unique solution of the ideal multiscale problem). Under Assumptions (Al)- 
(A3) and (B1 )-(B3), the ideal multiscale problem ([7]) has a unique solution. In particular, 
we have 


7(l + a l fi) 1 


< inf sup 

q&Qn vevjpi 


b(v,q) 

lkll Q ll v llv 


i.e. inf-sup stability independent of h and H. 



Proof. We let K™ h = {v e Vfff : V • v = 0}. The coercivity of a(-, •) on K™ h follows 
immediately from its coercivity on Kh since K™ h C Kh- The operator Id — Gh is stable 
in V with constant 1 + a” 1 /?, since V ■ = 0 and 

JK ; /> V II//■*($ i) < a _1 a(G h v, G h v) 

= a _1 a(v, G h \) 

< a 1 f 3 ||v|| i 2 (n) l|G h v|| 

L 2 (Q) 


for all v e V. Combining these results with the inf-sup stability of &(•, •) on Vh and Qh, 
we get 

b(v,q) 


7 < inf sup 


i£Qh v£ v h lkll Q ||v|| 


in 


< (1 + a 1 /3 ) inf sup 
= (1 + a _ 1 /3) inf sup 


(V ■ (Id — Gh)v,q) 
Q^Qh VGVh WqWq II (Id - G h )v\\ v 
(V • v, q) 


( 8 ) 


*QH^Jq llo l|v|| 


Qll V llV 


i.e. &(•, •) is inf-sup stable with constant y(l + independent of H and h. We note 

that K™ h = {v e Vjjf h : b(v,q H ) = 0 Wq G Qh}, since V ■ v <G Qh (see Remark [Ij). 
Finally, we apply Lemma [ 2 ] with V = Vfff, Q = Qh, /C = Kff h and constants a — a, 
(3 = (3 and 7 = 7(1 + a -1 /?)^ 1 . □ 


3.3 Error estimate for ideal problem 

In this section, we show that the flux solution of the ideal multiscale problem above 
converges in the energy norm with linear order in H to the reference solution. This 
convergence is independent of the variations of A, i.e. we do not have any pre-asymptotic 
effects from the multiscale features. 

Lemma 9 (Error estimate for ideal solution). Under Assumptions (A1)-(A3) and (Bl)- 
(B3), let u h solve Q and u 1 Jf i h solve (J7]). then 

IK - u“ HI < p^CaC^HWf - P H f\\v m 

where C p ^ and Cg are independent of h and H. 

Proof. Parametrizing the solutions u h(f) and u f} s h (f) by the data /, we use the triangle 
inequality to obtain 

|||w(/)-u^(/)||| 

< \\\n h (f)-n h (P H f)\\\ + \\\n h (P H f) -u- (P H f)\\\ + \\\u% h (P H f) -u^(/)|||. 

The two last terms will be shown to equal zero. 

For the first term, we proceed in several steps. Let us define := u h(f) — u h(Pnf) = 
u h(f — Ph f), which is the flux solution for the data / — Puf ■ The corresponding pressure 
solution shall be denoted by ph- First, we observe 

IIIW ||| 2 = (/ - Ph f,Ph) = (/ - Pnf,Vh ~ PHPh) < \\f ~ PHf\\i?<n)\\Ph ~ PHPh\\tf{n)- 

(9) 
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In order to bound the term ||ph — PHPh\\L 2 {Q), we let 0 G be the weak solution to 

A0 = ph — PHph■ Then we have 

= (Ph ~ PHph, 0 — Ph4>) A C p ^H\\p h -? PHPh 11z, 2 (n) 101ifi(r*)- 

Defining w := V0 we get V-w = p h -P H ph and ||w|| L 2 (n) < C Ptd H\\p h -P H p h \\ L 2 {n) . Next, 
we use a pair of projection operators 11^ : V —> 14 and Ph : Q —> Qh that commute with 
respect to the divergence operator, allows for Th to be non quasi-uniform, and where II/, 
is L 2 -stable, i.e. P/, V • w = V-If^w and ||IVw||£ 2 (q) < CgllwH^m), with C'p independent 
of h. The existence of such operators is proved in m Exploiting this stability and the 
fact that ph — PHph = i4(V • w) (since Ph is a projection on Qh and ph — PhPh £ Qh), 
we obtain 

II Ph ~ PHPh\\ 2 L 2 { a) = 0 Ph ~ PHph,Ph) = (-Ph(V • w ),p h ) 

= (V ■ U h w,p h ) = -(A _ 1 Uft,n fc w) < |||u /l |||||A- 1 / 2 n h w|| L 2 (a) 

< / 5 1/2 C' fi |||u fe ||||| w || L 2 (n ) < (3 1 / 2 CfjC p ^H\\\u h \\\\\p h - P H Ph\\L*(n)- 

Combining this estimate with (|9]) yields 

IKlf < P 1/2 C n C p4 H\\f - P H f\\ L 2 {a] \\\u h \\\. 

For the second term, since the correctors are divergence free, we have V ■ u'pf^P^f) G 
Qh- This implies V ■ ? h {P H f ) = — Ph f, hence 

V • u™ h (P H f) - V ■ u h (P H f) = 0, 

i.e. u™ h (P H f) — u h {P H f) G K h . Now, from first the equations in (|2j) and (J7|) in combi¬ 
nation with the a(-, •)-orthogonality between Vpp\ and K { h , we get 

a(u h (P H f), v) = 0, v G 14, V ■ v = 0, and 

a(u % h (P H f),v) = 0, v G 0™, V • v = 0, and 

v ) = 0, V G 14 f , V ■ V = 0. 

Since 14 = V44 © Tjf, we obtain 

a(u h (P H f) - u ™ h (P H f),v) = 0, 

for all v G K h . Choosing v = u h (P H f) - u ™ h (P H f), we see that u h (P H f) = u^ h (P H f), 
thus the second term equals zero. 

To show that the third term is zero, it is sufficient to show that VLpf h (f — Ph f) = 0. The 
data f—Pnf is L 2 -orthogonal to the test space Qh and it enters the equation ([T]) only in an 
L 2 scalar product with test functions. Hence u™ h {f—P h f) = u 40( Ph (,/— Ph f )) = 0. 0 

4 Localized multiscale method 

The ideal corrector problems ([5]) are at least as expensive to solve as the original reference 
problem. Hence, we require to localize these problems to very small patches, without sac¬ 
rificing the good approximation properties. If we can achieve this, the corrector problems 
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can be solved with low computational costs and fully in parallel. In this section, we show 
that this is indeed possible. We prove that we can truncate the computational domain 
in the local corrector problems (pH) to a small environment of a coarse element T. This 
is possible, since the solutions of (|5j) decay with exponential rate outside the coarse ele¬ 
ment T. We obtain a new localized corrector operator which can be used analogously to 
the ideal corrector operator to construct a localized multiscale space. This localization 
reduces the computational effort for assembling the multiscale space significantly. 

In addition to the assumptions (Al)-(A3) and (B1)-(B3), we require additional as¬ 
sumptions on the computational domain and the mesh. More precisely we assume the 
following for the analysis. 

(A4) We consider d — 2 and a simply-connected domain hi C M 2 . 

(B4) The fine grid Th is quasi-uniform, i.e. the ratio between the maximum and the 
minimum diameter of a grid element is bounded by a generic constant. 

We note that assumption (A4) is crucial for our proof. Assumption (B4) on the other 
hand could be dropped with a more careful analysis. In this case the estimates (and in 
particular the decay) will depend on the inverse of the minimum mesh size of the fine 
grid in a patch U(T). For simplicity of the presentation, we do not elaborate this case 
and restrict ourselves to quasi-uniform meshes, i.e. to (B4). Note that even though we fix 
d — 2, we keep the general notation d to illustrate how the results are influenced by the 
dimension. The localized method can be formulated analogously for d — 3. 

In order to localize the detail space K[ , we use admissible patches. We call this 
restriction to patches localization. For each T G Th we pick a connected patch U(T) 
consisting of coarse grid elements and containing T. More precisely, for positive G N we 
define fc-coarse-layer patches iteratively in the following way. For all T G Th (which are 
assumed to be closed sets), we define the element patch C4(T) in the coarse mesh Th by 


U 0 (T) := T, 

I'klT) (J{V € Th : T n Uk-i(T) 0} k= 1,2,.... 


( 10 ) 


See Figure [T] for an illustration of patches. For a given patch U(T), we define the 
restriction of V } [ to U ( T ) by 

I4 f ([/(T)) := {w G Vl : w = 0 in II \ U(T)}. 

Accordingly, we also define 

K[{U(T)) := {w G V*(U(T)) : V ■ w = 0}. 


Using this localized space, we define the localized corrector operators. Localized quantities 
are indexed by the patch layer size k. 

Definition 10 (Localized corrector operators). For each T G Th and k > 1 layers, we 
define a localized element corrector operator G kk : V —» K f h (Uk(T)): 


a(Gj fc v,w) = a T (v,w) 


( 11 ) 


for all w G K { h {Uk{T)). Further, we define the localized global corrector operator Gh,k ■ = 
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(a) One-coarse-layer patch, k = 1. 


(b) Two-coarse-layer patch, k = 2. 


Figure 1: Illustration of fc-coarse-layer patches. Dark gray subdomain is T. Light gray 
subdomain is C4(T). 


The localized corrector operators are again well-defined by the Lax-Milgram theo¬ 
rem, exploiting that a(-, •) is a weighted L 2 -scalar product. Note that the definition of 
Kh(Uk(T )) implies Neumann boundary conditions on the localized corrector problems 
(11). We define a localized multiscale function space by 


V£S fc :=(Id -G h , k )(V H ) 


and state the localized multiscale problem as follows. 


Definition 11 (Localized multiscale problem). The localized multiscale problem reads: 
find G V'lfili and pu 6 Qh, such that 


a ( u H S ’h> v h) + b(v h ,p H ) = 0, 


( 12 ) 


for all v h G Vfif k and q H G Q 


H- 


Definitions 10 and [IT] constitute the proposed multiscale method. Next, we show that 
the above stated problem is well-posed. 

Lemma 12 (Unique solution of localized multiscale problem). Under Assumptions (A1 )- 
(AS) and (Bl)-(BS), the localized multiscale problem (12) has a unique solution for all 
k, h and H. 

Proof. We use similar arguments as in Lemma |8j The basic difference is that we need to 
show stability for the localized corrector operator G^k- We start with the stability of the 
localized element corrector operators. Here we have for arbitrary v G V 


\ G h,k^\ 


= G h,k v ) = « T ( V > G h,k v ) < 


\ G h.,k^\ 


(13) 
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Now, we can prove Instability of the localized global operator. We get 


II^Wfc v llL 2 (n) — 


= a 


1 

- 2 °" 


Y G h,k v < «" 

TeT H L 2 (f2) 

1 Y a ( G h,k v , G h,k v 

TGT H T'cU k (T) 

1 E E 

TGT H T'<zU k (T) 


'« ( Y 

\T£Th T'gTh ) 


GLviir+ 


< fc v 


< oT l C p k d J2 |l|G^v||| 2 f a 'ey £ 


T < a f3C p k ||v|| i2(n) , 


T£T h 


T&Th 


where C p is a constant only depending on the shape regularity constant p of the coarse 
mesh. Similar to ([8]) we derive inf-sup stability with 

b(v,q) 


7 < inf sup 

q&Qn w&V H 


qIMIv 


< (1 + a~ lG f3 lG Cp G k d/2 ) inf sup 


q£Q H ve y”f 


!Qll v llV 


Observe that the inf-sup stability constant 7° := 7(1 + a l ^ 2 f3 l ^ 2 C G2 k d / 2 ) 1 depends on 
k this time. □ 

The inf-sup stability constant 7° depends on k due to overlapping patches. We come 


back to another estimate of the inf-sup stability constant in Section 4.3 after proving the 
decay of the correctors. 

It is important to note that in the localized case we do not have orthogonality between 
and K f h as in the ideal case (cf. equation Q). This orthogonality was crucial in the 
error estimate for the ideal method presented in Lemma [9j In the localized case, we rely 
on the exponential decay of the localized element correctors, which justifies localization 
to patches. 


4.1 Error estimate for localized problem 

In this section we state the main result of this paper, which is an a priori error estimate 
in the energy norm between the reference solution and the localized multiscale approx¬ 
imation. We first present a logarithmic stability result for the nodal Raviart-Thomas 
interpolation operator II h for fine scale functions and then state a lemma on the ex¬ 
ponential decay of the correctors. Then the main theorem follows. The proof of the 


exponential decay is contained in Section 4.2 The notation a < b stands for a < Cb with 


some constant C that might depend on d, a, (3 and coarse and fine mesh regularity 
constants, but not on the mesh sizes h and H. In particular it does not depend on the 
possibly rapid oscillations in A. 

We recall a well known stability result for the nodal Raviart-Thomas interpolation 
operator. 
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Lemma 13 (Logarithmic stability of the nodal interpolation operator for divergence free 
functions). Assume (B1 )~(BJ i ). For any given element T e Th there exists a constant C 
that only depends on the regularity of T and the quasi-uniformity of Th, such that 

ll n HV /t ||22 (T) < C , A(Lf/h)“||v/ ! ||^2( T ), 
with X(H/h) := (1 + log(H/ h)) 1 / 2 for all v h e 14 with V ■ = 0. 

A proof for this can be found in [HU, Lemma 4.1]. This result holds for both d — 2 
and 3. 


Remark 14. There exist unconditionally L 2 -stable Clement-type interpolation operators 
for which we could define X(H/h ) := 1 for all h and H instead, see 0 CZ3 EH [If- In 
particular, the operators introduced in are projections and were used as a technical 

tool in the proof of Lemma [5] above. However, these operators are hard to implement in 
practice and hence are not used in the proposed numerical method. 


Lemma 15 (Exponential decay of correctors). Under Assumptions (Al)-(Af) and (B1 )- 
(Bf), there exists a generic constant 0 < 6 < 1 depending on the contrast (3/a, but not 
on h or H such that for all positive k G N: 


( G 'I v - G 'L V ) 

T£T h 


< k d x{H/h) 2 e 2k /^ H/h) 

T&Th 



(14) 


for all v e V. 


Proof. The lemma is a direct consequence of Lemma 21 in Section 4.2 


□ 


Now, combining the error estimate for the ideal multiscale method in Lemma [9] and 
Lemma [15] we get the following a priori error estimate of the localized multiscale method. 


Theorem 16 (Error estimate for localized multiscale solution). Under Assumptions 
(Al)-(Aj) and (Bl)-(Bf), for a positive k E N, let u h solve (J2]) and solve (12), 
then 

HWf-PHfWmm + k^MH/h^/^WfW^a,, ( 15 ) 


ms,fc 
Uh — U Hh 


< 


for some 0 < 9 < 1 depending on the contrast fd/cx, but not on k, h and H. 


Before stating the proof, we discuss the role and choice of k. The second term in 
the error estimate (15) is an effect of the localization. This term can be made small by 
choosing large values of k, i.e. large patch sizes. A natural question is how to choose k to 
make the second term of order H to some power. 

We write A = A (H/h) for convenience. Let k = 2d~ 1 \og(6)X~ 1 k = — CgX~ 1 k, where 
Cq = —2 dr 1 log(6 l ) > 0 is a constant independent of H and h. We are interested in the 
asymptotic behavior, so we consider H 1. Setting the second term in (15) equal to 
HWfh^n) yields 

ke h = -C e X-^ d - 1 H 2/d , 


that is k = W(—CqX 4 / d 1 H 2 ^ d ), where W is the Lambert W-function. In terms of 
the number of layers k, we get k = — Cf 1 XW(—CeX~ A t d ~ l H 2 t d ). This equation has two 
solutions for sufficiently small H. Since we require k > 1, we pick the branch W < —1. 
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Another, more practical option is to choose k = R\\og(l/H) for some constant R. 
Then the expression k d ^ 2 X9 k ^ x will be asymptotically (as H —> 0) dominated by the power 
jj-Riog8 Choosing R sufficiently large yields arbitrary order of accuracy of the term. The 
fine mesh size h is often fixed and we can choose 

k = (1 + |log r .(i/)|) 1/2 log s (l/iJ) (16) 


for some bases r and s of the two logarithms. 

Remark 17. If Clement-type interpolation operators are used, we have X = 1 indepen¬ 


dent of H/h. Choosing k = C\og(l/H) makes the second term in (15) proportional to 


log(l /H) d / 2 H clo & e . For an appropriate C we can make the first term in (15) dominate 
the error estimate. 


Proof of Theorem \U\ Let uff£ := ((Id - G h>k ) o U H )uff h e Vfff K , then uffff - ufff is 
divergence free. Hence, by Galerkin orthogonality we have 




- ms,A: 


ms ,k • 


/ msi ms,Ac\ / ms,/c 

a {Uft - n H f , u h - n H h ) = a( u h - u H f , u h - u 


ms ,k\ _ 


ms ,k 


'ms,fe\ 
H,h ) 


and obtain 


ms,A; 

U h - U H,h 


< 


~ msi 
- U H,h 


< u h - U^J + 


..ms -ms.fc 
u H,h u H,h 


The first term can be bounded by ^^C^CpyHWf — Hh/Hl 2 ^) by Lemma [9j Regarding 
the second term, using [3~4 . Lemma 4.1] and stability of the ideal multiscale solution, we 
get 

£ |||Grn„u-||| 2 < £ lllnsu^lll 2 = |||n„uS* t ||| 2 <A(H/A) 2 ||/||i I(n) 

TgT h TgT h 

and can combine this with Lemma [15] to get 

ms -ms.fc 
u H,h u H,h 


= |||(G»,k-G k )n„u^||| 

£ (Gh - 


T&Th 


1/2 


< 


k d/2 X(H/h)9 k/xiH/h) HI GIXIhvlh) 

\T£Th 

<k d / 2 x{H/h) 2 e k ^ H ^\\f\\ L 2 m . 


□ 


4.2 Proof of exponential decay of correctors 


This section consists of four lemmas, Lemma 18 21, of which the last one is the main 


result. The two first lemmas are auxiliary and are motivated by steps in the proofs of 
the latter two. Before starting, we need to set some notation and introduce some tools. 
We use the notation = {/:/€ H l {pj) V compact subsets cu C M. d }. Note 

that we will use the letter K to denote arbitrary triangles of the coarse mesh 7 h- The 
first lemma says that every divergence free function w in H( div,fl) is the divergence of a 
skew-symmetric matrix. 
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Lemma 18. Let Q be a simply connected domain with Lipschitz boundary and let w G 
//(div, 12) with V • w = 0 in LI. Then there exists a skew-symmetric matrix if G 
[<’ c 2 (M d )] dxd with Vifij G [L 2 (M d )] d and J n 4> = 0 such that 

w = V-V> 12 and ||V^|| L 2 (w) < ||w|| L 2 (a;) for a; C 12. (17) 

Here, the divergence of if is defined along the rows. 

Note that the above lemma is the only instance, where we require the restriction d = 2. 
Even though the existence of the skew-symmetric matrix is also available for d — 3, we 
could not prove localized estimates of the type || V^qlln 2 ^) fS || w I|l 2 (w)- 


Proof. The result is a combination of well-known results. First, we extend the divergence- 
free vector field w G H(div, Q) to a divergence-free vector field w G H( div,M d ). In 
particular we have w G [Z/ 2 (M d )] d and w = w in Q. Note that the extension of w to 
will be typically not zero outside of 12. The existence of such an extension operator was 
proved in [S3J Proposition 3.8]. It is well known that there exists a skew-symmetric matrix 
if G [IT 1 Q C 2 (M d )] dxd with Vifij G [L 2 (M d )] d , such that w = V • if (see [23, Lemma 2.3]). 
The matrix is only unique up to a constant, so we fix the constant by J n if = 0 (which 
gives us a Poincare inequality). The inequality || V'0q||L 2 (aj) || w || (for uj C M d ) can 

be seen as follows for d = 2 . Obviously, if i — j we obtain Vifa = Tifj 3 = 0 and estimate 
01 is trivial. If i ^ j, we obtain by using the skew-symmetry 

II W IIl 2 (u;) = IIV • V’ll^H = ll^l^ll + d 2 lfl2\\ 2 L 2( ljJ ) + 11^1^21 + d 2 lf22 ||l 2 (u; ) 

= + ||<9l ifzi 11^2^) = + 11^1^121| 2 2^) 

= l|vVh 2 ||i2 H = ||v^2i||i2 H , 


i.e. 


we obtain even equality in estimate (17). 


□ 


We also require suitable cut-off functions that are central for the proof. For T G 77 
and positive k G N, we let the function i]T,k G P\{Th) (globally continuous and piecewise 
linear w.r.t. Th) be defined as 


VT,k(x) = 0 for x G £4_i(T), 
VT,k( x ) = 1 for x G 12 \ U k (T). 


(18) 


We start with the following lemma, which enables us to approximate truncated func¬ 


tions from Kf. 


Lemma 19. Let w h G Kf and let if G [W / ioc(^)] dX<i vn ^ 1 w ft, = V • if denote the cor¬ 
responding skew-symmetric matrix as in Lemma 18. Let furthermore ifx '■= \K\f K if 
denote the average on K G Th and let ifn G [L 2 (n)] dxd denote the corresponding piece- 
wise constant matrix with Vtff(^) = ifx for x G K. The broken divergence-operator 
Vh • is given by Vh • v := V ■ v \k for K G Th- The function r] T )fc G P\ (Th) is a 
given cut-off function as defined in (18) for k > 0. Then, we have that the function 
w h '■= lift (V • (r]T,k'f’)) ~ (n# o lift) (V ■ (r]T,k'f’)) G K f h fulfills the following estimate for 
any K G Th-' 


V ■ ( r]T,kf 0 — V h ■ (v t^h) ~ Wft|| L2(x) 

< f X(H/h)\\w h \\ L 2 {K) K C U k (T) \ U k ~i(T) 
~ ] 0 otherwise. 
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Obviously we also have supp(w/,) C0 \ Uk-i(T). 

Proof. First, we observe that the skew-symmetric matrix if must be a polynomial of 
maximum degree 2 on each fine grid element. We use this in the following without 
mentioning. 

We fix the element T G Th and fceN and denote 77 := r)T,k■ Furthermore, we define 
for I< G T„ 


c K := \K\ 1 [ rj and 0# := \K\ 1 f 0. 

J K J K 

We define w h := IT (V ■ (ijif)) — (11# o II#) (V ■ (ijif)) and observe that w h G K[ and 
w h = w h on tt\Uk(T). The property n#(w/0 = 0 is clear. The property V-w h = 0 follows 
from the fact that ?70 is still skew symmetric and that V • (II# o II#) (-) = (P# o P h ) (V ■). 
Since 0# and c# are constant on K we have 

n h (v • (c#0#)) = V • {c K if K ) = 0 on K. (19) 

Furthermore, since II#(v#) = v# for all v# G V# and since V ■ (r/0#) G V# we also have 

(n# o n A )(v ■ (#a 0 ) = n,, (v ■ (##)) o n K - ( 20 ) 

Finally, we also have on K , 


(n# o 110 ) (v • (c# 0 )) — c#(n# o n ft )(v • 0) — c#n#(w 0 — 0. (21) 


Combining (19), (20), and (21) we obtain for every K G 7# 


II(n# o n h ) (v • ( 770 )) - n h (v ■ (rnl> K ))\\v(K) 

= || (n# o n#) (v ■ (770) - v • (c# 0 ) - v ■ (770#) + v ■ (c# 0 #)) || L 2(#) 

= II (n# o n#) (v • ((77 — c #)(0 — 0#))) ||l 2 (a:)- (22) 


Now, we consider the quantity we want to estimate. For any K G 7#, 


V ■ ( 770 ) - V# ■ (770#) - w 0 | L 2 ( #) 

< || V • (77(0 - 0 #)) - II,, (V • (77(0 - 0 #))) || L 2 ( #) 

+ ||n# (V • (77(0 - 0#))) - n h (V • ( 770 )) + (n# o n ft ) (v ■ c#)) || L 2 ( # } 


= II v • (77(0 - 0#)) - n /t (v • (77(0 - 0#))) || L 2 ( #) 

+1|(n# o u h ) (v • (770)) - n,(v • ( 770 #))|| L 2 ( #) 

® IIv ■ ((77 - c #)(0 - 0#)) - n /t (v ■ ((77 - c# ) (0 - 0#))) 11 #2 (#) 

+ II (n# ° iih) (v • ((77 - c#) (0 - 0#))) ||l2 ( #) 

< X(H/h)\\V ■ ((77 - c #)(0 - 0#))|| L 2(#). ( 23 ) 


I11 the last step we used Lemma 13 , the property that II#V • ((77 — c#)(0 — 0#)) is 


divergence free and the fact that IT is locally L 2 -stable when applied to functions of 
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small fixed polynomial degree, i.e. for fixed t ETh and r G N there exists a constant C(r) 
that only depends on r and the shape regularity of t such that 


||n h (v)|| L 2 (t) < C(r) ||v|| £ a (t) for all v G [P r (t)] d . 


Continuing from (23) we obtain 


V- {{V ~ ck)^ - M\\h{K) 

£ II {V ~ c k )V • ijjf L 2 {K) + \\{ip- MVr]\\ 2 L 2 {K ) 

< ^llVrylll.^HV^Il!,^) 
f |l|w||i a(K) KcU^fU^T) 

~ 1 0 otherwise. 


(24) 


Note that we used the properties of r) to obtain the Lipschitz bound ||r/ — ck\\l°°(k) 
H\\Vri\\ L oo( K ) < 1 and that V?/ has no support outside C4(T) \ C4_i(T). We also used the 


Poincare inequality for // 
yields the sought result. 


ck which has a zero average on K. Combining (23) and (24) 

8 


We continue with a lemma showing the exponential decay of solutions to problems of 
the form in ([5]). 

Lemma 20. Now, let w 2 G K { h be the solution of 



■ V/, = F T (v h ) 


for all v h <E K[ 


(25) 


where F T G (N[)' is such that F T (v h ) = 0 for all w h G \ T). Then, there exists 

a generic constant 0 < 0 < 1 (depending on the contrast /3/a) such that for all positive 
k G N: 


w 


\n\u k (T) 


c Qk/\(H/h) I 


W 


T I 


In' 


(26) 


Proof. The proof exploits similar arguments as in j30|. L et us fix k G N. We denote again 
rj := r]T,k £ P 1 (Th) (as in (18)). We apply Lemma 19 to w 2 G K { h . The corresponding 
skew symmetric matrix shall again be denoted by if = -0(w 2 ) and we define 

w T : = n ft (v • (#)) - (n H o u h ) (v ■ (#)). 

We obtain that V ■ {ryi/f) — V h • (/IV h) — w 2 is zero outside Uk(T ) \ Uk-iiT) and 

IIV • (rifj) - V H ■ (vi’H) ~ w t || L 2 (t/ fc (T)\t/ fc _ 1 (T)) ^ \(H/h)\ |w T 11 l 2 {u k (T)\u k ^ (t))■ (27) 
First observe that 


A _ 1 w r ■ w 2 = / A _ 1 w 2 ■ w 2 ' = F t (w 2 ) = 0 


1 T 


’a \tr fc _i(T) 


(28) 


and 


r/w 2 = r/V • if — V • (r/VO — ybV/ 7 . 


(29) 
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With that we have 


A 1 w 7 • w 7 < 


>n\u k (T) 


'0\C/ fe _i(T) 


A~ 7 w 7 • (r/w 7 ) 


J29l 


ml 


A 1 w 7 • (V • (rjip) — -0V?/) 


'n\[4_i(T) 


ln\u k -i(T) 


'0\J7 fc i(T) 


A 1 w r ■ (V ■ — ifr'Vr] — w r ) 

A _1 w r • (V • ( 77 VO - • (iji/jh) ~ w T ) 


-V- 

=:I 


A 7 w T • (Vh ■ (#h) - V’Vry). 


' ^\U k —i(T) 


For I we use (27) to obtain 


I < X(H/h)\\\w 


— : 11 


T||| 2 

E/ fc (T)\t7 fc _i(T) 


and for II we obtain 

II = 

< 


A 7 w 7 ■ ((i> H - i/j)Vt}) 


'fiWk-dT) 

E 


w 


liC 


H\\'Vr]\\ L °a( K ' ) \\'V'ip\\L2{K) 


< 

nsj 


KgTh 

K<zU k {T)\U k ^(T) 

III T||| 2 
W 111 


I U k (T)\U k —i(T )' 

Now, denote by L := CX(H/h), and we get 


w 


T||| 2 


\n\u k (T) 


< L w 


T||| 2 


\U k {T)\U k ^(T) 


< L ( w 


.T || | 2 


I n\E4_!(T) 


w 


T||| 2 


I n\u k (T) 


where C is independent of T, k and A, but can depend on the contrast. We obtain 


.TUI 2 


,,I W HI Cl\U k (T) - ( X + L ■ , 

A recursive application of this inequality and w 7 


_lw lllw Tm2 


IO\t/ 0 (T) 


n Wfc-hT)' 

< lllw T | 


in 


yields 


w 


.Till 2 


\n\u k (T) 


< e 


-log(l+L-i)fc||| w T||| 2 ^ c -log(l+C?- 1 )fe/A(rr// ! ,)||| w T||| 2 


in — e 


w 


In’ 


where we used Bernoulli’s inequality and that 0 < L 1 < C 1 in the last step. The choice 
6 := (1 + G -1 ) -1 proves the lemma. □ 

The following lemma is the main result of this subsection. It can be directly applied 
to the localized corrector problems ( PH with F T (y h ) = a T (v,v h ), Gj lk v = w™ and 
= w 7 for any v e V. 
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Lemma 21. Let the setting of Lemma\20 \ hold true and let additionally w 1,k G K f h (U k (T )) 
denote the solution of 


'U k {T ) 


A 1 w T>k ■ v h = F T fv h ) for all v h G K f h (U k (T)). 


(30) 


Then, there exists a generic constant 0 < 6 < 1 (depending on the contrast) such that for 
all positive k G Kb 


^ ( w r _ w T ’ fc ) 

T&Th 


< 

r^j 


k d \(H/h) 2 e 2k/X{H/h) 

T&Th 


w 


12 
In' 


(31) 


Proof. Let r]T,k be defined according to (18) and denote z := ^ Tgrjf (w r — w i - fe ) G K { h . 
We obtain 

lll z lllo = 5 ^ (A -1 (w T — w T,fe ), (1 — r)T,k+i) z ) + (A“ 1 (w T «- w T,fe ), 77T,jfc+i z ). 


TsTh 

The first term is estimated by 

I < |||w T - w T ’ fe |||J|z(l - 7^r,fc+ 1 )||| I7 (T) < 


= :II 


|w T - w T ’ fc | 


InlH Hm fc+ i(T)- 


For the second term we have z G hence there exists again a skew-symmetric matrix 
if = if( z) with the properties as in Lemma 18 with 

VT,k+ 1 Z = ?7r,fc+iV • if = V ■ (r^fc+iVO - 


We define z := II/ l (V ■ (i]r,k+iff)) ~ (11# o IT#) (V ■ (r)T,k+iff))- Using Lemma 19 and 
supp( 77 T,fc+i z )nsupp(w r,fc ) = 0 we get 


(A x (w T - w r ’ k ),i]T,k+iz) = (A L w 1 , ? 7 r,fc+iz) 


IT 


} 28 l 


ln\u k (T) 


ln\u k (T) 


A w ■ (V ■ (r] T)k+l if) - ifVr)r tk+ 1 - z) 

A -1 (w r - w T,fc ) ■ (V • (rjT,k+i^) ~ -0V?/T,fc +i - z) 


Now proceed as in Lemma [20] to obtain 

II < X(H/h)\\\w T -w T ’ k \ 


lolll z IHc/ fc+ 1 (T)- 


Combining the estimates for I and II and applying Holder’s inequality finally yields, for 
k > l 

lll z llln X(H/h) |||w T -w T ' fe .~ m 


lnlN z, Hli/fc+ 1 (T) 


T&Th 


< k^\{H/h) | |||w 7 — w 

\T&T h 


(32) 


T,k 111 2 

In 


in- 
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It remains to bound |||w T — w T,fc |||^. In order to do this, we use Galerkiu orthogonality 
for the local problems, which gives us 

III T Tk\ l|2 ^ . r III T ~TA’II|2 

w — w ’ llL < mf w — w ’ L. 

w T ’ k &K{{U k (T)) m 


Again, we use Lemma 20 to show 


I T 

W — W 


T,fc||| 2 < 02k/\{H/h)^ 


W 


T||l 2 

In- 


Combining (32) and (33) proves the lemma. 


(33) 

□ 


4.3 Inf-sup stability revisited 

The decay results can be used to prove another inf-sup stability constant 7 ^ in addition to 


7 ° from Lemma 12 for the bilinear form &(■, •) with the localized multiscale space. Using 


Lemma 21. we obtain 


II G h ,kV - G h v 11^2(0) - 


- G h~ v ) 


TgTh 


T 2 (f2) 


< 


y ||Gjv||| 1(0) 
T&Th 

<k d MH/h) 2 0 2l ^ !h ^Ml Hnr 


We get the following stability 


||G\fcv|| L 2 (q) < || G h>k v — G h v\\ L 2 (ty + ||G'ftv|| 

< {k d / 2 \{H/h)9 k /* H ^ + l)||v|| i2(n) . 


12 


we obtain an inf-sup stability constant 7 ^. : = 


Using the same technique as in Lemma 

7(2 + k^XiH/h^WV)- 1 . 

For the nodal Raviart-Thomas interpolation operator II//, A (H/h) depends on h and 
H, and we cannot obtain a uniform bound on the constant for this estimate either. 


However, for L 2 -stable Clement-type interpolation operators (discussed in Remark 14), 
we have \{H/h ) = 1, independently of h and H. If using such an interpolator in place of 
n//, the inf-sup stability constant 7 ^ can be bounded from below by a positive constant 
independent of h and H , since k d ^ 2 6 k is bounded from above with respect to k. 


5 Numerical experiments 


Four numerical experiments are presented in this section. Their purpose is to show that 


the error estimate for the localized multiscale method presented in Theorem 16 is valid 
and useful for determining the patch sizes and that the method is competitive. 

A brief overview of the implementation of the method follows. The two dimensional 
Raviart-Thomas finite element is used. For all free degrees of freedom e (interior edges), 
the localized global corrector Gh^ e f° r the corresponding basis function <l> e is computed 
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according to equation 0 - The additional constraints on the test and trial functions to 
be in the kernel to the coarse Raviart-Thomas projection operator are implemented using 
Lagrange multipliers (in addition to those already there due to the mixed formulation). 
The corrector problems are cheap since they are solved only on small patches. This can 
be done in parallel over all basis functions. Finally, problem (12) is solved. Regarding the 
linear system arising here, we compare it with the linear system arising from a standard 
Raviart-Thomas discretization (using Vh for the flux) of the mixed formulation on the 
corase mesh: 

K B t ^ _ {O' 

B 0 


for matrices K and B and a vector b. The difference with the multiscale method is that 
matrix corresponding to the bilinear form a(-, ■) is computed using the low dimensional 
modified localized multiscale basis — Gh^E^E spanning V™ h k . Since the correctors 
are divergence free, K is replaced by a different matrix K in the system above, whereas 
B and b are left intact. 

In all numerical experiments below, the diffusion matrix is diagonal with identical 
diagonal elements, A(x) = A(x)I, with I being the identity matrix, for a scalar-valued 
function A. 


5.1 Investigation of error from localization 

In this experiment, we investigate how the error in energy norm of the localized multiscale 
solution is affected by the localization to patches of the correctors. The error due to 
localization is bounded by the second term in the estimate in Theorem 16 This term will 
be the focus of this experiment. 

The computational domain is the unit square = [0, l] 2 and the source function is 
given by 

f 1 if x e [0,1/4] 2 , 
f{x) = 1-1 if x 6 [3/4, l] 2 , 

I 0 otherwise. 


We consider three different diffusion coefficients A: 


1. Constant: A{x) = 1 in the whole domain. 

2. Noise: A{x) is piecewise constant on a 2' x 2 7 uniform rectangular grid. In each grid 
cell, the value of A is equal to a realization of exp(10cn), where a; is a cell-specific 
standard uniformly distributed variable. 

3. Channels: A(x) is as is shown in Figure [ 2 } It is piecewise constant on a 2' x 2' 
uniform rectangular grid. The coefficient A{x) — 1 for x in black cells and A{x) = 
exp(10) for x in white cells. 

Figure [3] shows the mesh used in the experiment. Both fine and coarse meshes are con¬ 
structed as shown in the figure. A reference solution u/j was computed with the standard 
Raviart-Thomas spaces 1 4 and Qh with h = 2 -8 . Solutions uj/Zf to the localized multi¬ 
scale problem were computed using H = 2~ 2 , 2~ 3 ,..., 2 -6 . The patch size k was chosen 
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Figure 2: Coefficient A defined on a 2 7 x 2 7 grid of Q. 


(a) Coarsest mesh, h = 1 . (b) One refinement, h = 1 / 2 . 

Figure 3: Family of triangulations of the unit square. 




as 

k = C(l + log 2 (H/h)) 1 / 2 log 2 (H- 1 ) 

rounded to the nearest integer with C = 0.25 and C = 0.5. The relative error (using the 
reference solution in place of the exact solution) in energy norm, i.e. u/, — u//)^ /|||uft||| 
was computed. See Figure [4] for the resulting convergence of this error with respect to 
H for the two values of C. Note that since / e Qh for all examples, the first term in 
(15) vanishes. The error is hence bounded by k d ^ 2 \{H/h) 2 0 k ^ x ^ H ^ h2> \\f\\L2^p j ^ which allows 
for a careful investigation of the influence of k, H and h. A reference line proportional 
to H 2 is plotted for guidance. We can see that we achieve convergence for both choices 
of C. However, since k is rounded to an integer, the convergence plots have a staggered 
appearance. This example shows that the error due to localization can be kept small and 
decreases with H. The plots also show the relative error in energy norm for the standard 
Raviart-Thomas discretization on the coarse mesh. It is evident that the localized mul¬ 
tiscale space has good approximation properties since it permits convergence while the 
standard space of the same dimension does not. 


5.2 Investigation of instability 


In this experiment we show how singularity-like features can appear in the solution, prob¬ 
ably as a result of high contrast in combination with the L 2 -instability of the nodal 
Raviart-Thomas interpolator. 

Again, we consider the unit square = [0, l] 2 . The diffusion coefficient A is chosen 
according to Figure [5j In other words, A is defined as 


A{x) 


exp( 10 ) if X 2 < 1/2 or x G [| — 2“ 5 , | + 2“ 5 ] x [§, \ + 2" 5 ], 
1 otherwise. 
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(a) Diffusion coefficient is constant. 



(b) Diffusion coefficient is noisy. 



(c) Diffusion coefficient has channel structures. 


Figure 4: Convergence plots for localization error experiments. Relative error in energy 
norm for three choices of A, for different values of the constant C determining the patch 
size. The number adjacent to a point is the actual value of k for the specific simulation 
corresponding to that point. 
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Figure 5: Coefficient A defined on a 2 5 x 2 5 grid of 12 where A(x) = 1 for x in black cells 
and A{x) = exp(10) for x in white cells. 


The source function is chosen as 


fix) = 


-1 if x 2 < 1/2, 
1 otherwise. 


This particular choice of A and / yields a localized multiscale solution with a clear 
singularity-like feature at x — (cc 1 , 0 : 2 ) = (1/2,1/2) in the localized multiscale solution. 

We use the family of triangulations presented in Figure [3] and fix H = 1/4 so that / 
is resolved on the coarse scale. Then / e Qh and all error stems from localization (see 
Theorem 16). We let the resolution h of the hire space be h = 2 -5 , 2 -6 ,..., 2~ 9 . Choosing 
k = 2, we compute the localized multiscale solution and reference solution u h for 

the given values of h. 

we expect to have 


From the error estimate in Theorem 16 


ms ,k 

u h - 


<k d / 2 x(H/h) 2 e k ^ H ^\\ f \\ L 2 m 

oc log (h _1 ) as h —* 0. 


The energy norm of the error is plotted in Figure [6} We can see that for this particular 
problem and range of h, the error increases with h and with the rate log(/r^ 1 ) as predicted 
by the error estimate. However, the error estimate seems not to be sharp for this particular 
example. Figure [7] shows the reference and multiscale flux solutions. The magnitude of the 
reference solution is in the range [0,3], while the multiscale solution has a spike reaching 
magnitude 30 at x = (1/2,1/2). Interesting to note is that the singularities vanish for the 
ideal multiscale method, i.e. without localization, see Lemma [9] 


5.3 Convergence in an L-shaped domain 

Next, we consider an L-shaped domain with noisy diffusion coefficient A (case 2. in Sec¬ 
tion 5.1) and with / ^ Qh- In this experiment, we show that the localization error 
investigated in the previous section can be dominated by errors from projecting /. 

We use the domain fl = [0, l] 2 \ [1/2,1] x [0,1/2] and the triangulation presented in 
Figure [8j Both fine and coarse meshes are constructed as shown in the figure. Further, 
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log 2 h 1 

Figure 6: Divergence of the energy norm of the localization error of a particular multiscale 
solution as h decreases. 



(a) Reference solution, h = 2 9 . (b) Multiscale solution, h = 2 9 , J7 = 2 2 , 

k = 2. 


Figure 7: Magnitude of flux at the centroid of the triangles. 


we choose the source function as 

( 1/2 + X 1 -X 2 if X 2 < 1/2, 

f( x ) = \ -(I/ 2 + xi - x 2 ) if x\ > 1/2, 

[ 0 otherwise. 

Note that / ^ Q# and ||/ — Ph/Wl 2 ^) ^ H. A reference solution was computed with 
the standard Raviart-Thomas spaces 14 and Qh with h = 2~ 8 . Solutions to the 

localized multiscale problem were computed using H = 2 -2 , 2 -3 ,..., 2 -6 . The patch size 
k was chosen as 

k = C(l + log 2 (H/h)) 1 / 2 log 2 (H- 1 ) 

rounded to the nearest integer, with C = 0.25 and C = 0.5. The relative error in energy 
norm was recorded for the solutions corresponding to the values of H. The resulting 
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(a) Coarsest mesh, h = 1/2. 


(b) One refinement, h = 1/4. 


Figure 8: Family of triangulations of the L-shaped domain. 


convergence plot can be found in Figure [9] We expect the first term in the error estimate, 


msi 

U h - U H,h 


< 


H\\f 


-Pff/|U 2 (n) 


(34) 


to be of order H 2 . From the convergence plots we can see that C = 0.25 is not sufficient 
to make the localization error of at least order H 2 , however, C = 0.5 is. 



“ C = 0.25 
•— ♦ C = 0.5 
- - oc H 2 


Figure 9: Convergence plot for experiment with L-shaped domain. Shows relative error 
in energy norm for two values of C and a series of values of H. The number adjacent to 
a point is the actual value of k for the specific simulation corresponding to that point. 


5.4 Comparison with MsFEM 

We compare the proposed method with the results obtained using the Multiscale Finite 
Element Method (MsFEM) based approach in [5], The domain is 12 = [0,1.2] x [0,2.2] 
and the permeability coefficient A is given in a uniform rectangular grid of size 60 x 220 
by the 85th permeability layer in model 2 of SPE10 [12]. 

The method proposed in [5j is based on a foie and a coarse mesh with quadrilateral 
elements. The fine mesh is uniform 60 x 220, i.e. aligned with the permeability data, 
and the coarse mesh is 6 x 22, so that each coarse element is subdivided into 10 x 10 
fine elements. The implementation of the method proposed in this work uses triangular 
meshes, which is why we divide each of the rectangular elements into two triangular 
elements by a diagonal line drawn from the upper left corner to the lower right corner. As 
coarse mesh, we use a similar triangular mesh that is constructed from a 6 x 22 rectangular 
mesh such that the fine mesh is a conforming refinement of the coarse mesh. 

The (quasi-singular) source data / is equal to 1 in the lower left and —1 in the upper 
right fine quadrilateral element. Note that such / is a discretization of point sources 
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that model production wells. In particular, the source terms on the continuous level are 
mathematically described by Dirac delta functions. Hence, for h —)• 0, we only have 
/ G W~ m,2 {Vt) for rn > |, opposed to / G L 2 (Q) as is required for our analysis. To 
account for this difference, we follow j28| and compute the localized source corrections 
F^ 1 f G V h l(Uc(T)) on t'-coarse-layer patches for T G Th, 


( 77 ’^\£ p 

a ( F h /> 


+ b(yl F^f) + b(F^f, q{) = -(/, g*) T , 


for all vl G I /£(Ug(T)) and G Q f h (Ue.(T )), where Q { h (Ue.[T )) is the restriction of Q { h 
to U(_(T ), analogous to the definition of V£(Ui(T)). (The pressure solution F^’ e f is not 
needed for correcting the flux and is discarded after its use as Lagrange multiplier). Since 
/ is non-zero only for the two triangles T\ and T 2 in the lower left and upper right corners, 
only two such corrector problems need to be solved. The total localized source correction 


is 


F hf = F f 


Ti,l 
h 


f + F h 


■T 2 ,l 


/e 


vl 


The localized corrector problems 0 are unaffected by the source correction. The 
right hand side of the localized multiscale problem (12) is appended with the localized 


source corrections and instead reads: find such that 


a ( u H S ,h’^ v h) + b{v h ,p H ) + b(u ™£'\q H ) = -(/, to) - a(F^f,v h ). 


ms.fci 


Using a value of t = 0 will be referred to as an ad-hoc source correction, since we do not 
expect to have any decay of the correction already within the source triangle itself. The 
source corrected solution is u^ fe,£ + F^f. 

We emphasize that the need for source correctors for singular source terms is not 
an exclusive drawback for our approach, but it is a common necessity shared by all 
comparable multiscale methods in this setting. In particular they are also used for the 
MsFEM-based approach in |5] that we use for our comparative study. 

The proposed localized multiscale method was used to solve for the flux in the de¬ 
scribed problem for three corrector patch sizes: k — 1,2, and 3. Three variants of source 
correction were used: i) without source correction, i.e. ii) with ad-hoc source cor¬ 
rection, i.e. + F[f for l = 0 (without interpolation constraint), and iii) with source 

correction, i.e. + F^f for £ = k : k + 1, 00 . A reference solution was computed 

on the fine mesh. Table [l] shows the relative energy norm and L 2 -norm of the difference 
between the localized multiscale solution and the reference solution for the different val¬ 
ues of k and l. The corresponding L 2 -norm of the error for the MsFEM method with 
oversampling HEO-OS proposed in [S] is also presented in the table. Note that HEO-OS 
is based on a discretization with roughly 33% less degrees of freedom than the proposed 
method, since it uses quadrilaterals instead of triangles (however, since this holds for both 
the fine and the coarse mesh, the relative change in the amount of degrees of freedom with 
respect to the reference solution is the same). The flux solutions are plotted in Figure |Toj 

The results show that the proposed method even without error correction compares 
favorably with the homogenization based approach. Ad-hoc error correction gives small 
errors for this problem in both norms. For source correction with patch size d = k, 
instabilities similar to that studied in Section A2 cause the error to increase. However, 
letting i = k + 1 is enough to get errors that compare favorably with [5]. 
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Tabic 1: Relative error in energy norm and L 2 -norm for the SPE10-85 problem. 


Method 

k 

£ 

Energy norm 

L 2 -norm 

Proposed method 

1 

— 

0.7863 

0.4069 

without source 

2 

— 

0.7856 

0.3369 

correction 

3 

— 

0.7855 

0.3325 

Proposed method 

1 

0 

0.1541 

0.2700 

with ad-hoc source 

2 

0 

0.1515 

0.1467 

correction (£ = 0) 

3 

0 

0.1537 

0.1379 


1 

1 

0.1090 

0.8292 


1 

2 

0.0459 

0.2703 

Proposed method 

1 

oo 

0.0350 

0.2504 

2 

2 

0.0549 

0.7453 

with source correction 

2 

3 

0.0185 

0.0517 

(£ — k, k + 1 , oo) 

2 

oo 

0.0150 

0.0490 


3 

3 

0.0080 

0.0178 


3 

4 

0.0051 

0.0424 


3 

oo 

0.0041 

0.0088 

HEO-OS 0 

— 

— 

— 

0.3492 
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(d)k = l,£ = l (e)k = l,£ = 2 (f)ife = 2,^ = 3 

Figure 10: Flux solutions for the SPE10-85 problem. Figure (a) shows the reference flux 
solution and (b-f) show the multiscale flux solutions for k — 1 and 2, and different source 
corrections (£ — — 1 means no source correction and £ = 0 means ad-hoc error correction). 
The color maps to the magnitude of the flux at the midpoint of the triangular elements. 
The colors map from 10 -5 (white) to 10~ 2 (black) and is saturated at white and black for 
lower and higher values, respectively. 
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